################################################################
## Simulation studies
################################################################

########################
## generated fake data
########################
set.seed(12345678)
N <- 300
df <- 3
epsilon1 <- rnorm(N) # norm
epsilon2 <- rt(N,df) # t with df 3
epsilon3 <- rchisq(N,df) # chisq with df 3
epsilon4 <- rald(N,0.3) # ald with tau 0.3
epsilon5 <- 2/3*rnorm(100) + 1/3*rnorm(100,sd = 1/10) # mixture of norm: kurtotic
epsilon6 <- 1/2*rnorm(100,-1,2/3) + 1/2*rnorm(100,1,2/3) # bimodal
epsilon7 <- 1/5*rnorm(100,-22/25,1) + 1/5*rnorm(100,-49/125,3/2) + 3/5*rnorm(100,29/250,5/9) # skewed

x <- runif(N,0,5)
qtau <- 0.3
beta0 <- 2
beta1 <- -1

pr <- 1 - pnorm(-beta0 - beta1*x)
y1 <- rbinom(N, 1, pr)
pr <- 1 - pt(-beta0 - beta1*x,3)
y2 <- rbinom(N, 1, pr)
pr <- 1 - pchisq(-beta0 - beta1*x,3)
y3 <- rbinom(N, 1, pr)
pr <- 1 - pald(-beta0 - beta1*x,0.3)
y4 <- rbinom(N,1,pr)
pr <- 1 - ( 2/3*pnorm(-beta0-beta1*x) + 1/3*pnorm(-beta0-beta1*x,sd=1/10) )
y5 <- rbinom(N,1,pr)
pr <- 1 - ( 1/2*pnorm(-beta0-beta1*x,mean=-1,sd=2/3) + 1/2*pnorm(-beta0-beta1*x,mean=1,sd=2/3) )
y6 <- rbinom(N,1,pr)
pr <- 1 - ( 1/5*pnorm(-beta0-beta1*x,mean=-22/25,sd=1) + 1/5*pnorm(-beta0-beta1*x,mean=-49/125,sd=1/10) + 3/5*pnorm(-beta0-beta1*x,mean=49/250,sd=5/9) )
y7 <- rbinom(N,1,pr)

xs <- data.frame(x = seq(0,5,length.out = N))

pr1 <- 1 - pnorm(-beta0 - beta1* (xs$x))
pr2 <- 1 - pt(-beta0 - beta1*(xs$x),3)
pr3 <- 1 - pchisq(-beta0 - beta1*(xs$x),3)
pr4 <- 1 - pald(-beta0 - beta1*xs$x,0.3)
pr5 <- 1 - ( 2/3*pnorm(-beta0-beta1*xs$x) + 1/3*pnorm(-beta0-beta1*xs$x,sd=1/10) )
pr6 <- 1 - ( 1/2*pnorm(-beta0-beta1*xs$x,mean=-1,sd=2/3) + 1/2*pnorm(-beta0-beta1*xs$x,mean=1,sd=2/3) )
pr7 <- 1 - ( 1/5*pnorm(-beta0-beta1*xs$x,mean=-22/25,sd=1) + 1/5*pnorm(-beta0-beta1*xs$x,mean=-49/125,sd=1/10) + 3/5*pnorm(-beta0-beta1*xs$x,mean=49/250,sd=5/9) )



##########################################################################
# probit models
##########################################################################
probitout1 <- glm(y1 ~ x,family = binomial("probit"))
probitout2 <- glm(y2 ~ x,family = binomial("probit"))
probitout3 <- glm(y3 ~ x,family = binomial("probit"))
probitout4 <- glm(y4 ~ x,family = binomial("probit"))
probitout5 <- glm(y5 ~ x,family = binomial("probit"))
probitout6 <- glm(y6 ~ x,family = binomial("probit"))
probitout7 <- glm(y7 ~ x,family = binomial("probit"))


probit_predict1 <- ifelse(predict(probitout1,type="response")>0.5,1,0)
probit_predict2 <- ifelse(predict(probitout2,type="response")>0.5,1,0)
probit_predict3 <- ifelse(predict(probitout3,type="response")>0.5,1,0)
probit_predict4 <- ifelse(predict(probitout4,type="response")>0.5,1,0)
probit_predict5 <- ifelse(predict(probitout5,type="response")>0.5,1,0)
probit_predict6 <- ifelse(predict(probitout6,type="response")>0.5,1,0)
probit_predict7 <- ifelse(predict(probitout7,type="response")>0.5,1,0)
xs = data.frame(x = seq(0,5,length.out = N))

pdf("figures/figure1_appendix.pdf",width=7,height=7)
layout(matrix(c(1,2,3,4,5,6,7,8,8), 3, 3, byrow = TRUE))
pred_prob1 = predict(probitout1,xs,type="response",se.fit=T)
plot(xs$x,pred_prob1$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="Normal distribution")
lines(xs$x,pred_prob1$fit+ -1.96* pred_prob1$se.fit,lty="dashed")
lines(xs$x,pred_prob1$fit+ 1.96* pred_prob1$se.fit,lty="dashed")
lines(xs$x,pr1,col="red")


pred_prob2 = predict(probitout2,xs,type="response",se.fit=T)
plot(xs$x,pred_prob2$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="T distribution (df=3)")
lines(xs$x,pred_prob2$fit+ -1.96* pred_prob2$se.fit,lty="dashed")
lines(xs$x,pred_prob2$fit+ 1.96* pred_prob2$se.fit,lty="dashed")
lines(xs$x,pr2,col="red")


pred_prob3 = predict(probitout3,xs,type="response",se.fit=T)
plot(xs$x,pred_prob3$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="Chi-square distribution (df=3)")
lines(xs$x,pred_prob3$fit+ -1.96* pred_prob3$se.fit,lty="dashed")
lines(xs$x,pred_prob3$fit+ 1.96* pred_prob3$se.fit,lty="dashed")
lines(xs$x,pr3,col="red")


pred_prob4 = predict(probitout4,xs,type="response",se.fit=T)
plot(xs$x,pred_prob4$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="ALD (tau=0.3)")
lines(xs$x,pred_prob4$fit+ -1.96* pred_prob4$se.fit,lty="dashed")
lines(xs$x,pred_prob4$fit+ 1.96* pred_prob4$se.fit,lty="dashed")
lines(xs$x,pr4,col="red")


pred_prob5 = predict(probitout5,xs,type="response",se.fit=T)
plot(xs$x,pred_prob5$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="Mixture of Normal distributions:\nKurtotic")
lines(xs$x,pred_prob5$fit+ -1.96* pred_prob5$se.fit,lty="dashed")
lines(xs$x,pred_prob5$fit+ 1.96* pred_prob5$se.fit,lty="dashed")
lines(xs$x,pr5,col="red")


pred_prob6 = predict(probitout6,xs,type="response",se.fit=T)
plot(xs$x,pred_prob6$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="Bimodal distribution")
lines(xs$x,pred_prob6$fit+ -1.96* pred_prob6$se.fit,lty="dashed")
lines(xs$x,pred_prob6$fit+ 1.96* pred_prob6$se.fit,lty="dashed")
lines(xs$x,pr6,col="red")


pred_prob7 = predict(probitout7,xs,type="response",se.fit=T)
plot(xs$x,pred_prob7$fit,type="l",xlab="X",ylab="Predicted Probability",bty="l",main="Skewed distribution")
lines(xs$x,pred_prob7$fit+ -1.96* pred_prob7$se.fit,lty="dashed")
lines(xs$x,pred_prob7$fit+ 1.96* pred_prob7$se.fit,lty="dashed")
lines(xs$x,pr7,col="red")

plot(NA, xlim=c(0,5),ylim=c(0,1),axes=F,
     xlab="",
     ylab="")
legend("top",legend=c("True predicted probability","Estimated predicted probability from Probit","95% Confidence Interval"),
       lty=c("solid","solid","dashed"),col=c("red","black","black"),cex=1.5,bty="n")

dev.off()

